//
clear;
//
Percorso = "/home/ciampa/Lavoro/Scilab/Lavoro/Sistemi Lineari/";
exec(Percorso + "/Metodi/EGPP.sci");
exec(Percorso + "/Metodi/SA.sci");
exec(Percorso + "/Metodi/SI.sci");
//
T_J = 0;
T_EG = 0;
NE = [];
//
N = 5000;
ripetizioni = 50;
           //  100:   TJ = 8.4d-3   TG = 4.2d-2   EJ = 1.6d-16   EG = 3.3d-17
           //  500:   TJ = 4.0d-2   TG = 2.7
           // 1000:   TJ = 7.9d-2   TG = 18.8
           // 2000:   TJ = 1.5d-1   troppa memoria necessaria per EGPP
           // 3000:   TJ = 2.3d-1   troppa memoria necessaria per EGPP
           // 4000:   troppa memoria per costruire A; senza costruire A:
           //         TJ = 3.1e-1
           // 5000:   senza costruire A:
           //         TJ = 3.9e-1
//
//A = 3*eye(N,N); A(2:N,1) = ones(N-1,1); A(1:N-1,N) = ones(N-1,1);
b = ones(N,1);
//
c = b/3;
//
x_esatta = ones(N,1)/6; x_esatta(1) = 1/4; x_esatta(N) = 1/4;
tic();
for j=1:ripetizioni,
    x = zeros(N,1);
    xn = x;
    for k = 1:35,
        s = -(x(1)+x(N))/3;
        for i = 2:N-1, xn(i) = s + c(i); end;
            xn(1) = c(1) - x(N)/3;
            xn(N) = c(N) - x(1)/3;
            x = xn;
            NE(k) = norm(x-x_esatta);
    end;
end;
T_J = toc()/ripetizioni;
//
for k = find(NE == 0), NE(k) = %eps; end; 
printf("\nerrore relativo J = %4.3e\n",norm(x - x_esatta)/norm(x_esatta));
RES = zeros(N,1);
RES(1) = 3*x(1) + x(N) - b(1);
for i = 2:N-1, RES(i) = 3*x(i) + x(1) + x(N) - b(i); end;
RES(N) = x(1) + 3*x(N) - b(N);
printf("\nmisura residuo = %4.3e\n",norm(RES)/norm(b));
printf("\ntempo per Jacobi (N = %d) = %4.3e\n",N, T_J);
//
scf(0);clf();plot([1:35]',log10(NE));xgrid();
//
tic();
//for j=1:ripetizioni,
    [S,D,P,info] = EGPP(A);
    [c,info] = SA(S,P*b,"disabilita");
    [x_EG, info] = SI(D,c,"disabilita");
//end;
T_EG = toc();
//
printf("\nerrore relativo EG = %4.3e\n",norm(x_EG - x_esatta)/norm(x_esatta));
printf("\ntempo per Gauss (N = %d)= %4.3e\n", N, T_EG);
//
